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Abstract 

We present a systematic study of the effectiveness of light quark mass reweighting. This method 
allows a single lattice QCD ensemble, generated with a specific value of the dynamical light quark 
mass, to be used to determine results for other, nearby light dynamical quark masses. We study 
two gauge field ensembles generated with 2+1 flavors of dynamical domain wall fermions with light 
quark masses mi = 0.02 (m n = 620 MeV) and mi = 0.01 (m n = 420 MeV). We reweight each 



ensemble to determine results which could be computed directly from the other and check the 
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consistency of the reweighted results with the direct results. The large difference between the 0.02 
and 0.01 light quark masses suggests that this is an aggressive application of reweighting as can 



be seen from fluctuations in the magnitude of the reweighting factor by four orders of magnitude. 



Never-the-less, a comparison of the reweighed topological charge, average plaquette, residual mass, 
pion mass, pion decay constant, and scalar correlator between these two ensembles shows agreement 
well described by the statistical errors. The issues of the effective number of configurations and 
finite sample size bias are discussed. An examination of the topological charge distribution implies 
that it is more favorable to reweight from heavier mass to lighter quark mass. 

PACS numbers: ll.15.Ha, 12.38.Gc 14.40.Be 
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I. INTRODUCTION 



Generating ensembles of gauge field configurations with dynamical quarks is usually the 
most expensive part of a lattice QCD calculation. Because of the renormalization of the 
quark mass, we are usually not able to identify in advance the input bare quark mass that will 
correspond to a target renormalized mass. Quark mass reweighting is a powerful technique 
that allows us to fine tune the sea quark masses to their physical or other desired values 
after an ensemble has been generated, avoiding the computationally expensive generation 
of new ensembles. As is explained below, a reweighted lattice quantity is computed by 
averaging the product of that quantity and a reweighting factor over the given ensemble. 
The reweighting factor is chosen to reproduce the effects of the change in the action that 
would result from changing from the simulated to the reweighted quark mass. 



Reweighting of the strange quark mass has 



Deen widely applied to accomplish such a fine 



tuning of the physical strange quark mass {l-jj or to study dependence on the dynamical 
strange quark mass Qj. Reweighting of the light quark mass has been used less frequently, 
presumably because the light quark masses used in most current lattice calculations are 
sufficiently heavy that reweighting methods may not be able to bridge the gap between 



the simulated and physical masses. However, Aoki et al. [5| have applied light quark mass 
reweighting to obtain results at the physical light quark mass. In the future, we expect 
there will be an increasing use of mass reweighting to adjust the light quark mass to its 
physical value. Another important application is the use of light quark mass reweighting 
to avoid difficulties caused by a small quark mass such as long autocorrelation times or a 
loss of numerical stability at small mass. For example in the case Wilson fermions, one 
intentionally performs a simulation at larger light quark mass and then reweights to the 



physical mass 



7J. This can also potentially reduce the computational cost. 



As an important and widely used technique in lattice QCD, quark mass reweighting has 
been carefully studied and examined, for example, using Wilson [6| or overlap fermions. 
However, to better understand how mass reweighting performs, we believe that addition 
study of the resulting statistical fluctuations and a systematic comparison of a reweighted 
ensemble with one directly generated without reweighting are needed. As has been empha- 



sized in earlier work such as Refs. 



8|, a small overlap in configuration space between the 



simulated ensemble and the target ensemble can lead to large fluctuations in the reweighting 



factors and unreliable results because of the lack of sufficient statistics. There are two issues. 
First, for a fixed sample size N, the reweighting factor is determined from what may be a 
biased estimator. (We will discuss this in detail in section IV). If the reweighting factors 
fluctuate too much, this bias may lead to a large systematic error. Second, without a suf- 
ficiently large number of configurations with non- negligible weight, the error estimate may 
not be reliable, or even when reliable, the estimated may be too large for the result to be 
useful. 

In this paper, we provide a systematic check of the reweighting technique by using two 
ensembles generated with the Iwasaki gauge action and the Domain Wall Fermion (DWF) 
action for the quarks. In order to see the limitations of reweighting, we consider a case 
were the method is used to achieve a large change in fermion mass. Specifically we reweight 
the light input quark mass from 0.02 (m n = 620 MeV) to 0.01 (m n = 420 MeV) and then 
check the reweighted results for many important physical quantities against those computed 
directly on the dynamical mi = 0.01 ensemble. We also do this in the reverse direction, 
reweighting the mi = 0.01 ensemble to 0.02. With the ensemble size up to 1000 configura- 
tions (a total of 10,000 molecular dynamics time units), we find that even for such a large 
range of quark mass reweighting, the reweighted results and those obtained directly on the 
dynamically generated ensemble are consistent, demonstrating the success of the reweighting 
technique. 

We study the behavior of reweighting on both short and long distance quantities. The 
paper is organized as follows: We begin in Sec. Hi] with a description of the reweighting 
method, a detailed study of the reweighting factors and an analysis of the effective number 
of configuration that results after reweighting. We then discuss the reweighting of the 
topological charge and susceptibility in Sec. IIHI In Sec. IIVt we discuss the reweighting of a 
short distance physical quantity, the average plaquette, and the problem of reweighting bias. 
A second short distance quantity which shows significant dynamical quark mass dependence, 
the residual mass, is examined in Sec. [V] as well as two important physical observables: the 
pion mass and decay constant,. Finally the effect of reweighting on the scalar correlation 
function is presented in Sec. |\T]and in Sec. IVIII we further discuss and summarize our results. 
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II. LIGHT QUARK MASS REWEIGHTING METHOD 
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We follow closely the strategy presented in Ref [6|. Reweighting an ensemble of config- 
urations from the light sea quark mass mi to m 2 requires the evaluation of a reweighting 
factor w(U ; mi, m 2 ) for each configuration on which measurements are to be performed. We 
can see that the expectation value for any operator O on an ensemble with light sea quark 
mass m 2 can be calculated from an ensemble generated with sea quark mass m-i as follows. 
One begins with the usual expression for the expectation value of the operator O computed 
for the case of a sea quark mass m 2 : 

det{D\m 2 )D(m 2 )} y/det{D^(m s )D(m a )} 
det{Dt(l.0)D(1.0)} > /det{Dt(1.0)£)(l = 0)7' 

Here D(m) represents the DWF Dirac operator for a fermion of mass m, m s is the strange 
quark mass and the determinant factors with mass m — 1.0 are the standard Pauli-Villars 
contributions which appear in the DWF formulation. The quantity Z 2 appearing in the 
denominator is the partition function for the sea quark mass m 2 : simply the integral in 
Eq. ([1]) with the factor in the numerator representing the operator O removed. 

Next one multiplies and divides by factors of Z\ and det{D^(U, mi)D(U, mi)} obtaining 

(0} 2 (2) 

det{D\m 2 )D(m 2 )Y 



Z 2 Z\ 



o 



det{Dt(mi)D(mi)} 



Sg d^{D\m x )D[m x )} ^/det{D^(m s )D(m s )} 
9 det{£>t(1.0)D(1.0)} y/fet{Dt(1.0)D(1.0fi 

= (Ow(m ll m 2 )} 1 
{w{m 1 ,m 2 )) 1 

Here we have introduced the reweighting factor w(U ; mi, m 2 ) for each configuration U : 

det{D\U,m 2 )D{U,m 2 )} 

w(U ; mi, m 2 = , fnUTT r^rryj 4 

det{D'(U, mi)D(U, mi)} 

and performed a similar manipulation to write the ratio Z\jZ 2 as (w(mi,m 2 )) 1 . In the 
above formulae, we have used the notation D(U, m) to make explicit the dependence of the 
DWF Dirac operator on the gauge field variable U. The reweighting factor w(U;mi,m 2 ) 
can be determined from the stochastic average: 

f £)£ e -{e + £ , (mi)1'i3(m2) t - l D(m2)- :L D(mi)e-^g} e -^t^ 

w(U ; mi, m 2 ) = j (5) 

-{CtD(m. 1 )tD(m 2 )t-iD(m. 2 )- 1 D(m 1 )C^t 5 }^ ^ ^ 



where the average appearing in Eq. fl6]) is to be performed over the stochastic variable £ 
drawn from the random Gaussian distribution exp{— 

We have purposely arranged the product of the four D(U,m) operators in the exponent 
of Eq. ([6]) in a form that will minimize the cost of the needed Conjugate Gradient (CG) 
inversion. Let us call the exponent of that equation f(U ; mi, 1712; £) so that w(U; ml,m2) = 
( e -/(J/;mi,m2;0\ . Then /([/; mi, m 2 ; £) = rfrj - where 

rj = D(m 2 )- 1 J D(mi)£ = [ J D(m 2 ) t D^)]" 1 D{rw^ D(m x )£. (7) 

If m 2 is close to mi, then the solution for rj will usually converge faster than will the 
solution for a combination such as [D(m 2 )^i}(m 2 )] _1 -D(mi)£, which would have resulted 
from a different arrangement than that used above. 

It was pointed out in Ref j(| that the average of stochastic sources and the average of 
gauge fields can be performed in any order, so we can perform the average needed in Eq. (jf)]) 
using any number of random hits iV hit . By doing so, we introduce extra noise into our 
system, but no systematic error. To reduce the noise following the idea of determinant 
factorization proposed by Hasenbush , we calculate the reweighting factor by carrying 
out iV s tep steps. For each step we change the mass by Am = (m 2 — mi)/N step . Following this 
approach, the reweight factor for a given configuration is then calculated from the product 

iVstep 

w(U; ml, m2) — j [ w(U, mi + (i — l)Am, mi + iAm) (8) 

s=l 

iVstep r iVhit ^ 

_ J~J J e -/([/;mi+(i-l)Am,mi+iAm;{i, a ) I _ ^ 

8=1 I ^ hit i=l J 

where for each step, A^it hits are performed by summing over the iVhit Gaussian variables 
{£i,s}i<i<iv hit for each of the steps s, 1 < s < N stcp . 

Our goal is to calculate the reweighting factor with sufficient accuracy that the uncer- 
tainty introduced by the finite number of the random hits is smaller than the fluctuation 
of w(U;ml,m2) among the underlying gauge configurations. The total number of DWF 
Dirac inversions required by Eq. is N t = A^itA^tep. By increasing A" step , we reduce the 
fluctuation among the independent steps. Averaging over an additional Ahit hits gives a 
better estimate of the average. From our experiments with a fixed number of total inver- 
sions, we find that smaller errors can be achieved from more steps than more hits. These 
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results confirm similar conclusions reached by the RBC and UKQCD collaborations when 
reweighting in the strange quark mass. 

This observation can be verified by the following explicit calculation. Consider a single 
gauge configuration and evaluate the reweighting factor w as a product of the reweighting 
factors w s , 1 < s < N stcp corresponding to iV st0 p steps. In our stochastic evaluation the factor 
w s corresponding to the step s is given by a function w s (l;) of a set of random variables £ 
where a different stochastic vector £ will be used for each step. Define the quantity a s from 
the ratio of averages: 



= ^/N Lp (1Q) 

where the average is taken over the stochastic vector £ in the fixed gauge background. The 
factor 1/jVjj. has been extracted from the exponent on the right hand side because the 
exponent in Eq. ([6]) will be reduced in size by 1/N step when a fixed mass interval is divided 
into N stcp segments. Since the reweighting factor w is linear in each of the independent 
stochastic factors iy s (£ s ), it is a straight forward exercise to express the fluctuations in w in 
terms of the quantities a s . If w is obtained from a single hit: 

(5wf = ((w - (w)f) (11) 

= <> 2 > - H 2 (12) 

-/Vgtep -/Vstep 

= n k 2 > - n w a ( i3 ) 

8=1 8=1 
'Nstep 

T 2 /AT2 



{w) 2 f n e^/^-p - lj (14) 
(w) 2 (e a2 / N ^ - l) , (15) 



where in the final line we have made the simplifying assumption that cr s is independent of s. 
If this process is repeated N^ it times and the results averaged (iV^t hits) then the resulting 
error from the stochastic sampling becomes 



(Sw) = (w) x l^— , (16) 

A smaller error results if the iVhit hits are averaged at each step, as is indicated in Eq. 02]). 
For this case a series of steps similar to those in Eqs. (fTTT)-([15l) gives: 

1/2 




5w = (w)Ul + —(e"*»-l)\ -1} . (17) 



where we write the simpler expression that results if we continue to use uniform steps and 
assume that the corresponding quantities a s are independent of s. 

While the fluctuations given in Eq. ( JT7|) are smaller than those in Eq. ffl6|) . both ex- 
pressions are decreasing functions of N step for a fixed total number N t of Dirac inversions, 
making it favorable to use larger iV ste p instead of larger N^ t . When N step is sufficient large, 
both of the above formulae become approximately wa/y/Wf in this limit it is unimportant 
how N t is divided between steps and hits. A second advantage of introducing a number of 
Hasenbusch steps is that we can then reweight to many different intermediate mass values 
in a single reweighting calculation. A third advantage can be seen from Eq. ([7]): splitting 
the ratio m 2 /mi into more factors makes the mass difference appearing in each step smaller 
so that each inversion converges faster. Thus, in all our reweighting calculations, we use 
only one hit, but many intermediate mass steps. 

We use two ensembles generated with 2+1 flavors of domain wall fermions, the Iwasaki 
gauge action with (3 = 2.13, a 16 3 x 32 space-time volume and the extent L s = 16 in 
the fifth dimension.. The lattice spacing for this (3 value has been determined to be 



a -1 = 1.73 (3) GeV 10] . The configurations we use begin with trajectory 160 and extend to 
trajectory 3155 sampled every 5 trajectories for the m/ = 0.02 ensemble (600 configurations 
in total) and extend from trajectory 500 to 10490, sampled every 10 trajectories, for the 
mi = 0.01 ensemble (1000 configurations in total). We reweight the mi = 0.02 ensemble to 
mi = 0.01, and mi = 0.01 to mi = 0.02 through 80 steps performing only one hit for each 
configuration. The calculated reweighting factors are shown in Fig. [TJ 

Because of the large shift in mass that we are attempting to achieve with reweighting, 
the reweighting factors presented in Fig. [T]vary in size by a few orders of magnitude. Such 
large fluctuations in the reweighting factors may cause concern that such a substantial 
reweighting is doomed to fail. However, the large fluctuations of these factors toward lower 
values are not important. No mater how small they are, their only consequence is that 
the corresponding configurations make a negligible contribution in the reweighted ensemble. 
By simply counting the number of configurations that have reweighting factor larger than 
1, we see a variation in size of 100 and 200 for the 0.02 and 0.01 ensembles respectively. 
Thus, the reweighting may be more successful than Fig. [1] naively suggests. The efficiency 
of reweighting can be estimated in a more systematic way as follows. 

Suppose that we have a specific ensemble of N configurations. The reweighted value for 
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FIG. 1. Left: normalized reweighting factors computed from the mi = 0.02 configurations which 
can be used to reweight the light quark mass from m; = 0.02 — > 0.01. Right: normalized reweighting 
factors computed from the mi = 0.01 configurations which can be used to reweight from mi = 
0.01 — > 0.02. Here we plot the normalized reweighting factors vbi = Nw%/ '^2jWj which become 
unity in the case of equal weights. 

any physical quantity O based on Eq. (J3J) is 



(O w ) 



N 



(18) 



With the assumption that the reweighting factors Wi and the measurements O; are weakly 
correlated, we can derive an approximate expression for the statistical fluctuations of the 
reweighted average (O w ) N around the average value (((O w ) N )). Here the brackets ((...)) 
indicate the result that would be obtained when many ensembles of N configurations are 
generated and the results of each such calculation averaged. Of course, such an average will 
normally give the true value, i.e. the value in the limit of large N, except for quantities 
which are not simple averages where a bias may result if N is insufficiently large. If we 
keep only the dominant term we find (for a more detailed discussion and a derivation, see 
Appendix C of Ref. Q): 



(O w ) N - (((O w ) N )) 



50 



2 ^corr 



w- 



N w 



(19) 



where the quantities appearing on the right-hand side can be approximately calculated from 
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a single ensemble of N configurations according to the equations: 

N 



5 ° 2 = ^T,^-(°)n) 2 (20) 



i=i 



^=^i> ( 2i ) 

8=1 

— 1 - 

w2 =nY. w " ( 22 ) 
i=i 

^max 

r C o rr = ]T c(/)iy(/). (23) 

i = -trnax 

The quantities C(/) and W(l) appearing in Eq. f[2"31 are the autocorrelation functions for 
the measured quantity, O, and the the reweighting factors ty, and can be estimated from a 
single ensemble using the formulae: 

c ^-j^j2^ jo* { } 

i=l 

j=i 

for the case / > 0. When evaluating the autocorrelation time with a finite sample size N, we 
must impose a maximum length L max JV. Here, we use the smallest value of I at which 
C(l + 1) becomes negative. 

If we interpret the error given in Eq. ([19]) as the size of the fluctuations among our 



samples divided by the square root of an effective number of configurations, a/50 2 /N c g, we 
can determine iV e fT to be: 

N w 2 

N cS = (26) 

'corr W 

In the case of no reweighting(ifj = 1), this expression reverts to our normal expression for the 
effective number of configurations N/t cott . Even though different physical quantities have 
different autocorrelation times, the reweighted auto correlation times r corr in Eq. (123]) are 
usually very close to 1 for our data, because C(l) and W(l) are both small numbers except 
when 1=0. Therefore, the estimated effective number of configurations can be calculated from 
the reweighting factors alone. We find N c g = 48 for the mi = 0.02 ensemble (reweighting to 
0.01) and iV e fr = 63 for the m; = 0.01 ensemble (reweighting to 0.02). Clearly, the effective 
sample size decreases dramatically because of our large range of mass reweighting. However, 
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with a sample size of roughly 50 configurations, we can expect that the reweighting still 
works well for many interesting physical quantities. 

To justify our assumption that r corr ~ 1, we calculated the autocorrelation function C(l) 
for a few of quantities that we will study later. Figure [2] shows the autocorrelation function 
C(l) of the average plaquette(P), topological charge(Q) and susceptibility (Q 2 ) as well as 
the autocorrelation function W(l) defined in Eq. (j25p . The W(l) function converges to the 
constant value w 2 /w 2 for large separations. From these results the integrated autocorrelation 
times r corr defined by Eq. fl23|) for the unreweighted and reweighted ensembles are calculated 
and summarized in Table [TTJ We can see that even for the topological charge which has the 
longest autocorrelation time, the reweighted integrated autocorrelation time is less than 2. 
Since most quantities have significantly shorter autocorrelation time compared to that of the 
topological charge, we can expect that r — > 1, thus supporting our estimate of the number of 
effective configurations. As can be seen in Tab. HI1 the reweighted ensembles typically show 
a decreased autocorrelation time. This should be expected because reweighting effectively 
thins the initial ensemble. Note this decrease in r corr partially compensates for the decrease 
in iV e fj arising from the U 2 /u 2 factor inEq. (1261) . 




10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 

trajectory trajectory 

FIG. 2. The autocorrelation function C{1) defined in Eq. (I24p plotted versus I for the average pla- 
quette (P) , topological charge (Q) and Q 2 . The autocorrelation function W(l) for the reweighting 
factor defined in Eq. (|25|) is also shown. The left panel shows results calculated from the mi = 0.02 
ensemble while that on the right shows the results from the mi = 0.01 ensemble. 
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TABLE I. The integrated autocorrelation time defined in Eq. (|23p for both the unreweighted and 
reweighted data. Notice that the time is given in units of 5 trajectories for the mi = 0.02 ensemble 
and 10 trajectories for the mi = 0.01 ensemble. The columns are labeled by the light quark mass 
used to generate the original ensemble. 



unreweighted 


mi = 0.02 


mi = 0.01 


reweighted 


mi = 0.02 


mi = 0.01 


P 


1.8 


1.1 


p 


1.1 


1.0 


Q 


7.5 


3.8 


Qw 


1.8 


1.1 


Q 2 


3.8 


2.5 


Ql 


1.4 


1.1 



III. TOPOLOGICAL CHARGE AND SUSCEPTIBILITY 

The topological charge typically has the longest autocorrelation time among the usual 
observables and represents long distance physics. We will examine the reweighting effect 
on it first. We will later study short distance quantities such as the average plaquette and 
residual mass. 



As the first step in determining the topological charge we use Wilson flow 111, [12] devel- 
oped by Liischer to smooth or cool the gauge configurations. As a function of the Wilson 
flow time, the gauge fields obey the following differential equation and initial conditions: 

dUti ^ /i) = iglY,j- s S w {e^ sXa U t )\ s=0 T a U t {x^), U t (x,^)\ t=0 = U(x,^), (27) 

a 

where S W (U) = p- ^ p ReTr{l — U(p)} is the Wilson gauge action and {T a }i< a < 8 the eight 
hermitian generators of SU(3). When acting on a particular link matrix U(y, v) the matrix 
X a = T a if y = x and v — \i and is zero otherwise. We follow the Runge-Kutta integration 
scheme (see Appendix C of Ref [3] ). This realizes a continuous version of the stout smearing 



technique developed by C. Morningstar and M. Peardon [13(. We have used an integration 
time step 5t = 0.04, which gives an error smaller than 10~ 5 for t < 16 as seen by comparing 
the results to those obtained with St = 0.005 for dozens of configurations. All results given 
here use St = 0.04 and the integration error should be negligible for the topological charge 
calculation of interest. 

We then determine the topological charge using the "5Li" (5-loop-improved) algo- 
rithm [l4| on these "cooled" configurations, finding values of topological charge very close 
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to integers. Figure [3] shows a few typical evolution curves for the topological charge Q as a 
function of the flow time t. Notice that there are ambiguities for some of the configurations: 
the topological charge converges to one integer, stays there for a while and then moves away 
converging to a neighboring integer. The problem of defining the topological charge on the 
lattice and its ambiguities are not the focus of this paper. We will simply adopt the Wilson 
flow prescription and use the value of Q at t = 16 as our result for the topological charge. 




-6 1 1 1 1 1 1 1 1 1 

2 4 6 8 10 12 14 16 

t 



FIG. 3. Topological charge Q(t) as the function of the Wilison flow time t. We show results from 
those configurations numbered between 500 and 4500 separated by steps of 100 trajectories from 
the mi = 0.01 ensemble. 

Figures H] (a) and (b) show the evolution of the topological charge for the mi = 0.02 and 
mi = 0.01 ensemble respectively. There are apparent autocorrelations among configurations 
which we have also seen in the corresponding autocorrelation functions shown in Fig. [2j The 
lower two graphs of Fig. H] present the topological charge distributions and the corresponding 
distributions obtained after reweighting. The reweighted topological charge distribution is 
given by the probability function P W (Q) computed from the relation: 

P W {Q) = (28) 

where the Kronecker delta factor Sq.^q ensures that sum over %' only includes those config- 
urations having topological charge Q. 

The most important property of the topological charge distribution (given its central 
value of zero) is the width (Q 2 ) of the distribution, which is related to the topological 
susceptibility. The calculated values of (Q 2 ) with the average (Q) and the reweighted values 
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(Qw) and (Q w ) are summarized below the distributions shown in Fig. HJ The quoted errors 
are calculated from the fluctuations seen among blocks of 100 trajectories to remove any 
effect of autocorrelation. 

We can see very clearly how we have successfully reweighted one distribution to another. 
If we reweight the light sea quark mass from 0.02 to 0.01, the resulting mi — 0.01 ensemble 
reweighted from mi = 0.02 gives a distribution with second moment (Q 2 W ) = 5.9(1.4), in ex- 
cellent agreement with the 'true' value from our dynamically generated mi = 0.01 ensemble, 
which gives (Q 2 ) = 6.17(38). In fact, a detailed comparison of the two mi = 0.01 ensembles 
in Fig. H] (the reweighted histogram in the bottom row, left column and the directly com- 
puted histogram in the middle row, right column) shows very similar features. Comparing 
the tails of the distribution one clearly sees that the configurations with large topological 
charge are strongly suppressed by the reweighting factors: the reweighting really reflects the 
underlying physics by giving small weights to large topological charge configurations and 
large weights to small topological charge configurations. Similarly, the reweighting from the 
mi = 0.01 ensemble to that with mi = 0.02 is also successful but looks less impressive as can 
be seen by comparing the reweighted histogram in the bottom row, right column and the 
directly computed histogram in the middle row, left column. This should be expected as 
we try to reweight a narrower distribution to a wider one. The reweighting must give large 
weights to the configurations with large topological charge, which are rare in the mi = 0.01 
ensemble. Thus, in the reweighted distribution the far ends of the tails which are popu- 
lated in the directly computed ensemble are missing from the reweighted one: these bins are 
empty in our finite sample and it is not possible to create something by reweighting. This 
may suggest that it is more favorable to reweight from heavier mass to lighter mass. 

The spike at Q = —5 in the mi = 0.02 distribution reweighted from the mi = 0.01 
ensemble in Fig. H] (bottom row, right-hand column) may look troublesome. We find that 
its contribution comes to a large extent from a single configuration that has the largest 
reweighting factor in Fig. [TJ Therefore, this large spike has very large inherent uncertainly. 
To evaluate the quality of the data set, we (incorrectly) dropped this configuration and 
performed a complete analysis, repeating all the calculations in this paper. This experiment 
revealed no surprises, giving results consistent with those presented here which included all 
configurations. Nonetheless, the unphysical topological charge distribution for the ensemble 
reweighted to mi = 0.02 ensemble suggests that our calculation is not far from being sensitive 
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FIG. 4. Topological charge evolution (top row) for the mi = 0.02 ensemble, measured every fifth 
trajectory (left), and for the mi = 0.01 ensemble, measured at every tenth trajectory (right). The 
middle row shows the resulting topological charge distribution for mi = 0.02 (left) and m/ = 0.01 
(right). The bottom row shows the histograms after reweighting: mi = 0.02 — > 0.01 (left) and 
mi = 0.01 — > 0.02 (right) . Shown below the histograms are the results for (Q) and (Q 2 ) for the 
original and reweighted ensembles. 
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m, 

FIG. 5. Reweighted values for the average plaquette as the dynamical light quark mass is decreased 
from mi = 0.02 to mi = 0.01 using the mi = 0.02 ensemble (open circles). Also shown are results 
from the m; = 0.01 ensemble as mi is increased by reweighting from mi = 0.01 to mi = 0.02. The 
data from the mi = 0.02 ensemble are slightly offset in order to show the overlapped data points 
more clearly. While each reweighting sequence was performed in 80 steps, we plot only the result 
of each eighth step of 0.001. 

to uncontrolled statistical fluctuations. The reweighting is not ideal. It would be prudent 
either to have more configurations or to reweight less aggressively in the mass difference. 

IV. AVERAGE PLAQUETTE 

Now let us discuss how well reweighting works when applied to a representative short 
distance quantity, the average plaquette. The reweighting is performed in 80 steps, so we 
can examine all the intermediate reweighted results. We plot the reweighted results from 
each mass change of 0.001 for reweighting in both directions in Fig. |5j As we can see, the 
average plaquette is successfully reweighted from 0.587938(18) for the mi = 0.02 ensemble 
to the 0.588103(49) value predicted for mi = 0.01 as well as from the 0.588047(11) for 
the mi = 0.01 ensemble back to 0.587971(28) value predicted for m; = 0.02. The results, 
including all the intermediate mass values, agree within la. The quoted error is calculated 
using a block size of 50 trajectories to remove any effect of autocorrelation. We will also use 
a block size of 50 trajectories to calculate the error in our later calculation of m res , /„. and 
m n . At this point, we have shown that light-quark mass reweighting is very successful for 

15 




cr 
as 



0.58816 
0.58814 
0.58812 
0.5881 
0.58808 
0.58806 
0.58804 
0.58802 
0.588 
0.58798 




cr 

03 



0.58804 
0.58803 
0.58802 
0.58801 
0.588 
0.58799 
0.58798 
0.58797 
0.58796 
0.58795 
0.58794 




1000 



10 



100 
N size 



1000 



FIG. 6. Biased estimator of the reweighted average plaquette. For mi = 0.02 ensemble(left), we 
choose sample size to be N=10, 20, 30, 50, and 100. The corresponding number of measurements 
are 60, 30, 20, 12 and 6. For mi = 0.01 ensemble(right), we choose sample size to be N=10, 20, 
50, 100, and 200. The corresponding number of measurements are 100,50,20,10 and 5. The data 
is then fit with function f{x) = P w (l + c/N). 

both long distance and short distance physics. 

So far we have ignored the bias which arises from the finite sample size used in our 
stochastic estimate of the reweighted quantity (O w )n given by the ratio in Eq. ( TT8l) . Let us 
take the average plaquette as an example and see how large it is. Suppose the sample size is 
N. The reweighted result (P w )j\[ is biased in the sense that the expectation value {{{Pw)n)) 
will not agree with the true value P w = lhxiTv^oo (P w )n- It is straightforward to work out 
the difference, keeping only the lowest order terms: 

Eli w i p i 



«<^>AT» 



N N 



P,A i 



L^yy 

N 2 w 2 P 



WiWi(P w - Pi 




(29) 



(30) 



i=l 1=1 

We can see that the bias is proportional to 1/N, and the coefficient is related to the 
correlation of W{ and Pj. We might attempt to estimate the bias given in this equa- 
tion for our single sample by omitting the average over samples, ((•••)) and by examining 
the / = % term which might be expected to be dominant. This term can be written as 
— Eij=i w i w j( w i ~ w j)(Pi ~ Pj)/2N 3 w 3 . This expression hints that the bias is negative for 
a positive correlation between Wi and Pi, and vice versa. As should be expected, if there 
is no correlation between Pi and Wi, the coefficient vanishes. Taking all of our data and 



1(3 



dividing them into small samples, we calculate the expectation value of (P W ) N . The results 
are shown in Fig. [6] and suggest that the bias for the mi = 0.02 ensemble is negative and the 
bias for the mi = 0.01 ensemble is positive although these effects are at or below the one 
u level. These results also suggest that the bias is negligible for reasonably large ensemble 
size N > 20. 



V. PION MASS AND DECAY CONSTANT 



Let us next study the effect of reweighting on some additional, interesting physical quan- 
tities. In this section, we will examine the sea quark dependence of the residual mass, the 
pion mass, and the pion decay constant. In the following, for the mi = 0.02(0.01) ensem- 
ble, we fix the valence quark mass at 0.01(0.02) and gradually reweight the sea quark mass 
from 0.02(0.01) to 0.01(0.02), and study the change in each of these physical quantities. 
Once each measurements has been reweighted to mi = 0.01(0.02), it becomes unitary with 
equal valence and sea light quark masses and we compare the reweighted result with the 
unitary result directly calculated from the mi = 0.01(0.02) ensemble. Note we will express 
all quantities in lattice units unless other units are explicitly specified. 
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FIG. 7. Reweighted results for the residual mass. The blue triangular (red circle) points are the 
residual mass calculated with mi = 0.02(0.01) ensemble, using valence quark mass 0.01(0.02). The 
two box-shaped points give the results from a direct, unitary calculation with m va i = 0.02 on 
mi = 0.02 ensemble, and m v . A \ = 0.01 on mi = 0.01 ensemble. They are intentionally shifted by a 
small amount horizontally to avoid overlap with the reweighted results. 



The residual mass is an important quantity in the DWF formulation, representing a 
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contribution to the physical quark mass from the explicit chiral symmetry breaking which 
is present if the lattice extent in the fifth dimension is finite. (For this calculation the fifth 
dimensional extent has the value L s = 16.) The residual mass is an interesting quantity to 
study because it typically shows a strong dependence on the light sea quark mass and so 
should provide an interesting test of light sea quark mass reweighting. Our results obtained 



using the methods of Ref. [lOj for the residual mass from both the mi = 0.02 and mi = 
0.01 ensembles are shown in Fig. [7J This figure shows that the residual mass is positively 
correlated with the light sea quark mass (at fixed valence quark mass). We see that the 
ensemble with sea quark mass rrt; = 0.02 when reweighted to mi = 0.01 gives a result that 
is consistent with that obtained by direct calculation on the m ; = 0.01 ensemble, and vice 
versa. 

To obtain m w and /„-, we calculated three different correlation functions. Using a gauge- 
fixed pseudo-scalar wall source, we choose the sink to be either a gauge-fixed pseudo-scalar 
wall sink, a pseudo-scalar point sink or a point sink formed from the zeroth component of 
the axial current. These correlation functions are identified as Cq^^ OUVCC2 (t) which for our 
three cases becomes C]^? / (t), Cp^{t) and C^p(t). Here the label £ W means a gauge-fixed 
wall source or sink, 'L' means a point sink, 'P' means pseudo-scalar operator, 'A' means the 
zeroth component of the axial current operator. We perform a simultaneous fit to all three 
correlators at large times, using the functional form 

Co\o 2 (t) = N^e-^ + e ~ m ^ T ^} (31) 
This fit give us m n directly. From the amplitudes Nq*q : we construct /„. following Ref j^j], 



2 (N™y 



f*- Z A\\- — mww ( 32 ) 



m w Np P 



where is the renormalization constant of the local axial current A^, calculated in Ref. 15 ]. 
Since Za is to be evaluated in the chiral limit and thus does not depend on sea quark mass, 
we simply set Za = 1 for the present calculation. 

Figure El shows the calculated dependence of the pion mass on the light sea quark mass. 
In the left graph, we show the reweighting of the mi = 0.02 ensemble, using a fixed valence 
quark mass m va i = 0.01. The results from the intermediate reweighting steps are also 
shown. Similarly, the right graph shows the reweighting on the mi = 0.01 ensemble. These 
results demonstrate a positive dependence of the pion mass on the light sea quark mass. 
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FIG. 8. The pion mass computed on the mi = 0.02 ensemble (left panel) at a fixed valence quark 
mass m va _i = 0.01 as a function of the reweighted light sea quark mass with the results shown as 
open circles. The right panel shows similar results for the pion mass computed using the mi = 0.01 
ensemble for fixed m va i = 0.02. The triangular points in each panel show the result computed on 
the other ensemble with equal sea and valence quark masses: mj = 0.01 (left panel) and mi = 0.02 
(right panel). 
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FIG. 9. The reweighted pion decay constant plotted in the same style as used in Fig. [HJ 



The reweighted results show a clear trend, although with increasing errors, as we reweight 
the sea quark mass. In both cases the result obtained in the final step of the reweighting 
procedure agrees with the directly calculated result at the 1 a level. Figure. |9] shows similar 
results for /„., presented in the same way as the data for m n . in Fig. |HJ Similar conclusions 
can be drawn. 
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VI. SCALAR CORRELATOR 



Finally, let us look at something less successful but equally interesting. The properties 
of the iso-vector, scalar particle are known to be difficult to calculate because its mass is 
much larger than the mass of the pion, so large are statistics are needed. For the reweighted 
data, the effective number of configurations drops to less than one hundred, resulting in 
such a small ensemble size that it is not possible to determine the mass of this scalar state. 
However, by simply looking at the scalar correlator itself, we can see the dramatic effect of 
reweighting: the partially quenched violation of positivity disappears as the sea quark mass 
is reweighted to the unitary value. 

We calculate the scalar correlator using an identical, gauge-fixed wall source and sink so 
that the resulting correlation function should be positive. A well known artifact of partial 
quenching appears in the case where m va i < m sea : the correlator becomes negative for 
large source sink time separations when the lighter, negative norm irrj' intermediate state 



dominates 



161 ] . This behavior can be clearly seen in the left panel of Fig. [TU] where the 



valence quark mass is 0.01, and the sea quark mass is 0.02. The correlator becomes negative 
for t > 4. After reweighting the sea quark mass from 0.02 down to the unitary value of 0.01, 
the correlator data points increase and becomes positive as seen previously in Ref. In 
the left-hand panel one can also see the result from a unitary calculation performed directly 
with m va i = 0.01 on the mi = 0.01 ensemble. These results agree within the larger statistical 
errors of the reweighted data. In contrast, the partially quenched result with m va i = 0.02 
determined on the mi = 0.01 ensemble gives a positive value for the correlator as shown in 
the right panel of Fig. [10J The reweighting from m sea = 0.01 to 0.02 now has a less dramatic 
effect but does give a result in agreement with a direct calculation with m va i = 0.02 on the 
mi = 0.02 ensemble. Because of the poor statistics, the comparison between the reweighted 
and direct results is not accurate, especially at large t where the relative error increases 
dramatically and the results become unreliable. 



VII. SUMMARY AND DISCUSSION 

We have shown how reweighting can be used to vary the light sea quark mass while 
working with a single ensemble generated with a single sea quark mass. We first successfully 
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FIG. 10. Reweighting results for the scalar correlation function plotted as a function of time. The 
left panel shows this propagator computed with a valence quark mass of 0.01. The open circles 
give the non-unitary result from 0.02 ensemble, the crosses the same observable reweighted to 
the unitary, light sea quark mass of 0.01 and the triangles the result of a direct calculation on 
the mi = 0.01 ensemble. The right panel shows the corresponding reweighting on the mi = 0.01 
ensemble. 

reweighted the topological charge distribution from one of our light sea quark masses to 
the another. We found that it is better to reweight from a heavier mass to a lighter mass, 
because this corresponds to reweighting from a wider distribution to one that is more narrow. 
In this case, the poor statistics present in the tail of the distribution is suppressed instead 
of magnified by the reweighting factor. Next we studied the reweighting of the average 
plaquette, examining results obtained for a sequence of intermediate masses. In addition, 
we examined the bias introduced by our stochastic finite-sample reweighted measurement 
of the average plaquette. We found that the bias is negligible for fairly large ensemble with 
N > 20. After this study of both long distance and short distance quantities, we examined 
the behavior under reweighting of the residual mass, pion mass and pion decay constant. In 
each case, reweighting revealed how the given physical quantity varied with the sea quark 
mass when computed on an ensemble with a fixed quark ass. In each case, the final step of 
reweighting was chosen to duplicate the sea quark mass for which had been used to generate 
a second ensemble. This allowed a direct comparison between the reweighted and a directly 
generated result and in each case the two results agreed within errors. Finally, reweighting 
was shown to dramatically to correct the partial quenching artifacts seen in the calculation 
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of the iso-scalar, scalar correlator initially computed with m VSL \ = 0.01 on an mi = 0.02 
ensemble. 

The error on the quantities we have studied was generally increased 3-4 times by reweight- 
ing, in rough agreement with our estimate of the number N e g of effective configurations. 
These estimates suggested an error increase of a factor of ^600/48 = 3.5 and ^1000/62 = 4, 
if autocorrelation is neglected. The topological charge shows the longest autocorrelation 
times and the associated errors do not follow our estimated value of N c q. However, the 
relative error on Q 2 increases by a factor of 1.5 when reweighting the mi = 0.02 ensemble 
and by 3 for the mi = 0.01 ensemble. As explained in the discussion Sec. [Til this possibly 
smaller increase in error may result from there being fewer independent configurations in the 
original ensemble so that less information is lost when a portion of these correlated samples 
is assigned a small weight. 

It should be emphasized that the change in the quark mass from mi = 0.02 to 0.01 and 
the reverse is a large change to be accomplished by reweighting. This is certainly suggested 
by the large fluctuations of more than 4 orders of magnitude seen in the reweight factors. 
The resulting efficiency in the reweighting procedure, N e s/N, is less than 10%, causing the 
error after reweighting to increase by a factor of 4. However, such an aggressive use of 
reweighting may be appropriate for a test of the method, more clearly revealing its power 
and potential limitations. 

In real applications of reweighting, we may choose to be more conservative and reweight 
over a smaller range of mass than considered here to keep the efficiency high. Had we 
only reweighted by half of the mass difference, reweighting both ensembles to mi = 0.015, 
the effective number of configurations would increase dramatically from 48 to 308 for the 
mi = 002 ensemble, and from 63 to 386 for the mi = 0.01 ensemble. Depending on the 
physics being studied, we must ensure that iV c fj is large enough to give statistically useful 
results. For example, N e s < 100, is apparently insufficient for a study of the iso- vector, 
scalar particle mass. In conclusion, reweighting has been shown to work remarkably well, 
even for a large change in the light quark mass, with the resulting errors well described by 
an effective number of configurations given in Eq. [26j 
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